Linear Mixed Models

Single Random Effect

data(reeds, package="brinla")
summary(reeds)
 site     nitrogen   
 A:5   Min.   :2.35  
 B:5   1st Qu.:2.62  
 C:5   Median :3.06  
       Mean   :3.03  
       3rd Qu.:3.27  
       Max.   :3.93  
library(lme4)
mmod <- lmer(nitrogen ~ 1+(1|site), reeds)
summary(mmod)
Linear mixed model fit by REML ['lmerMod']
Formula: nitrogen ~ 1 + (1 | site)
   Data: reeds

REML criterion at convergence: 13

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.221 -0.754 -0.263  0.914  1.612 

Random effects:
 Groups   Name        Variance Std.Dev.
 site     (Intercept) 0.1872   0.433   
 Residual             0.0855   0.292   
Number of obs: 15, groups:  site, 3

Fixed effects:
            Estimate Std. Error t value
(Intercept)    3.029      0.261    11.6
library(INLA)
formula <- nitrogen ~ 1 + f(site, model="iid")
imod <- inla(formula,family="gaussian", data = reeds)
summary(imod)

Call:
"inla(formula = formula, family = \"gaussian\", data = reeds)"

Time used:
 Pre-processing    Running inla Post-processing           Total 
         0.9781          0.1516          0.2377          1.3674 

Fixed effects:
              mean     sd 0.025quant 0.5quant 0.975quant   mode   kld
(Intercept) 3.0293 0.1402     2.7587   3.0293     3.3004 3.0293 1e-04

Random effects:
Name      Model
 site   IID model 

Model hyperparameters:
                                         mean     sd 0.025quant 0.5quant 0.975quant   mode
Precision for the Gaussian observations 13.25  5.203      5.471    12.49      25.59 10.940
Precision for site                      20.35 25.120      1.933    12.81      84.76  5.093

Expected number of effective parameters(std dev): 1.243(0.5604)
Number of equivalent replicates : 12.07 

Marginal log-Likelihood:  -20.88 
invsqrt <- function(x) 1/sqrt(x)
sdt <- invsqrt(imod$summary.hyperpar[,-2])
row.names(sdt) <- c("SD of epsilson","SD of site")
sdt
                  mean 0.025quant 0.5quant 0.975quant    mode
SD of epsilson 0.27471    0.42753  0.28294    0.19767 0.30233
SD of site     0.22168    0.71916  0.27942    0.10862 0.44309
prec.site <- imod$marginals.hyperpar$"Precision for site"
prec.epsilon <- imod$marginals.hyperpar$"Precision for the Gaussian observations"
c(epsilon=inla.emarginal(invsqrt,prec.epsilon),site=inla.emarginal(invsqrt,prec.site))
epsilon    site 
0.29086 0.31354 
sigma.site <- inla.tmarginal(invsqrt, prec.site)
sigma.epsilon <- inla.tmarginal(invsqrt, prec.epsilon)
c(epsilon=inla.mmarginal(sigma.epsilon),site=inla.mmarginal(sigma.site))
epsilon    site 
0.26648 0.22175 
inla.contrib.sd(imod)$hyper
                                    mean       sd    2.5%   97.5%
sd for the Gaussian observations 0.28960 0.058561 0.19695 0.42688
sd for site                      0.31705 0.165926 0.11014 0.76824
sampvars <- 1/inla.hyperpar.sample(1000,imod)
sampicc <- sampvars[,2]/(rowSums(sampvars))
quantile(sampicc, c(0.025,0.5,0.975))
    2.5%      50%    97.5% 
0.086604 0.495696 0.898860 
library(brinla)
bri.hyperpar.summary(imod)
                                    mean      sd  q0.025    q0.5  q0.975    mode
SD for the Gaussian observations 0.29079 0.05833 0.19817 0.28283 0.42608 0.26647
SD for site                      0.31313 0.15740 0.10893 0.27909 0.71371 0.22176
alpha <- data.frame(imod$marginals.fixed[[1]])
library(ggplot2)
ggplot(alpha, aes(x,y)) + geom_line() + xlim(2,4) + geom_vline(xintercept = c(2.66, 3.40))+xlab("nitrogen")+ylab("density")

x <- seq(0,1,len=100)
d1 <- inla.dmarginal(x,sigma.site)
d2 <- inla.dmarginal(x,sigma.epsilon)
rdf <- data.frame(nitrogen=c(x,x),sigma=gl(2,100,labels=c("site","epsilon")),density=c(d1,d2))
ggplot(rdf, aes(x=nitrogen, y=density, linetype=sigma))+geom_line()

bri.hyperpar.plot(imod)

sdres <- sd(reeds$nitrogen)
pcprior <- list(prec = list(prior="pc.prec", param = c(3*sdres,0.01)))
formula <- nitrogen ~ f(site, model="iid", hyper = pcprior)
pmod <- inla(formula, family="gaussian", data=reeds)
pmod$summary.fixed
              mean      sd 0.025quant 0.5quant 0.975quant   mode        kld
(Intercept) 3.0293 0.30703     2.4026   3.0293     3.6569 3.0293 2.2118e-07
bri.hyperpar.summary(pmod)
                                    mean       sd  q0.025    q0.5  q0.975    mode
SD for the Gaussian observations 0.28651 0.055529 0.19794 0.27906 0.41494 0.26358
SD for site                      0.46599 0.219836 0.18247 0.41706 1.02812 0.33867
bri.hyperpar.plot(pmod)

(lambda <- 3*sdres/tan(pi*0.99/2))
[1] 0.022066
halfcauchy <- "expression:
              lambda = 0.022;
              precision = exp(log_precision);
              logdens = -1.5*log_precision-log(pi*lambda)-log(1+1/(precision*lambda^2));
              log_jacobian = log_precision;
              return(logdens+log_jacobian);"
hcprior <- list(prec = list(prior = halfcauchy))
formula <- nitrogen ~ f(site, model="iid", hyper = hcprior)
hmod <- inla(formula, family="gaussian", data=reeds)
bri.hyperpar.summary(hmod)
                                    mean       sd  q0.025    q0.5  q0.975    mode
SD for the Gaussian observations 0.28763 0.056198 0.19809 0.28006 0.41769 0.26436
SD for site                      0.42336 0.231601 0.14712 0.36566 1.03014 0.28080
reff <- pmod$marginals.random
x <- seq(-1.5,1.5,len=100)
d1 <- inla.dmarginal(x, reff$site[[1]])
d2 <- inla.dmarginal(x, reff$site[[2]])
d3 <- inla.dmarginal(x, reff$site[[3]])
rdf <- data.frame(nitrogen=x,density=c(d1,d2,d3),site=gl(3,100,labels=LETTERS[1:4]))
ggplot(rdf, aes(x=nitrogen, y=density, linetype=site))+geom_line()

bri.random.plot(pmod)

sdres <- sd(reeds$nitrogen)
pcprior <- list(prec = list(prior="pc.prec", param = c(3*sdres,0.01)))
formula <- nitrogen ~ f(site, model="iid", hyper = pcprior)
pmod <- inla(formula, family="gaussian", data=reeds, control.compute=list(config = TRUE))
psamp <- inla.posterior.sample(n=1000, pmod)
psamp[[1]]
$hyperpar
Precision for the Gaussian observations                      Precision for site 
                                10.2559                                  6.2126 

$latent
              sample1
Predictor:01  2.90317
Predictor:02  2.89826
Predictor:03  2.90155
Predictor:04  2.90375
Predictor:05  2.90598
Predictor:06  2.65167
Predictor:07  2.64810
Predictor:08  2.65079
Predictor:09  2.64952
Predictor:10  2.65175
Predictor:11  3.46508
Predictor:12  3.46871
Predictor:13  3.47142
Predictor:14  3.46862
Predictor:15  3.47009
site:A        0.21266
site:B       -0.03975
site:C        0.77777
(Intercept)   2.69090

$logdens
$logdens$hyperpar
[1] -5.4707

$logdens$latent
[1] 72.867

$logdens$joint
[1] 67.397
lvsamp <- t(sapply(psamp, function(x) x$latent))
colnames(lvsamp) <- row.names(psamp[[1]]$latent)
mean(lvsamp[,'site:C'] > lvsamp[,'site:A'])
[1] 0.992

Longitudinal Data

data(reading, package="brinla")
ggplot(reading, aes(agegrp, piat, group=id)) + geom_line()

library(lme4)
lmod <- lmer(piat ~ agegrp + (1|id), reading)
summary(lmod)
Linear mixed model fit by REML ['lmerMod']
Formula: piat ~ agegrp + (1 | id)
   Data: reading

REML criterion at convergence: 1868.7

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.946 -0.595 -0.028  0.507  2.868 

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 29.8     5.46    
 Residual             44.9     6.70    
Number of obs: 267, groups:  id, 89

Fixed effects:
            Estimate Std. Error t value
(Intercept)  -11.538      2.249   -5.13
agegrp         5.031      0.251   20.04

Correlation of Fixed Effects:
       (Intr)
agegrp -0.949
formula <- piat ~ agegrp + f(id, model="iid")
imod <- inla(formula, family="gaussian", data=reading)
imod$summary.fixed
                mean      sd 0.025quant 0.5quant 0.975quant     mode        kld
(Intercept) -11.5352 2.26057   -15.9781 -11.5353    -7.0964 -11.5352 8.0785e-13
agegrp        5.0306 0.25287     4.5335   5.0306     5.5272   5.0306 1.0001e-12
bri.hyperpar.summary(imod)
                                   mean      sd q0.025   q0.5 q0.975   mode
SD for the Gaussian observations 6.7096 0.35619 6.0457 6.6959 7.4437 6.6650
SD for id                        5.2988 0.62982 4.1389 5.2738 6.6074 5.2391
bri.hyperpar.plot(imod)

summary(imod$summary.random$id$mean)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -7.869  -2.934  -0.145   0.000   2.645  14.017 
bri.random.plot(imod)

reading$cagegrp <- reading$agegrp - 8.5
lmod <- lmer(piat ~ cagegrp + (cagegrp|id), reading)
summary(lmod)
Linear mixed model fit by REML ['lmerMod']
Formula: piat ~ cagegrp + (cagegrp | id)
   Data: reading

REML criterion at convergence: 1819.8

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.654 -0.541 -0.150  0.385  3.281 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 id       (Intercept) 35.72    5.98         
          cagegrp      4.49    2.12     0.83
 Residual             27.04    5.20         
Number of obs: 267, groups:  id, 89

Fixed effects:
            Estimate Std. Error t value
(Intercept)   31.225      0.709    44.0
cagegrp        5.031      0.297    16.9

Correlation of Fixed Effects:
        (Intr)
cagegrp 0.563 
nid <- length(levels(reading$id))
reading$numid <- as.numeric(reading$id)
reading$slopeid <- reading$numid + nid
formula <- piat ~ cagegrp + f(numid, model="iid2d", n = 2*nid) + f(slopeid, cagegrp, copy="numid")
imod <- inla(formula, family="gaussian", data=reading)
imod$summary.fixed
               mean      sd 0.025quant 0.5quant 0.975quant    mode        kld
(Intercept) 31.2241 0.69723    29.8520  31.2241    32.5947 31.2241 1.1108e-12
cagegrp      5.0305 0.29162     4.4567   5.0305     5.6036  5.0305 1.4099e-12
bri.hyperpar.summary(imod)
                                    mean       sd  q0.025    q0.5  q0.975    mode
SD for the Gaussian observations 5.52166 0.305437 4.94884 5.51143 6.14773 5.49058
SD for numid (component 1)       5.68547 0.557750 4.66364 5.65933 6.85172 5.61227
SD for numid (component 2)       1.90716 0.265745 1.44226 1.88667 2.48399 1.84545
Rho1:2 for numid                 0.94699 0.043666 0.82917 0.95948 0.98841 0.97515
postmean <- matrix(imod$summary.random$numid[,2],nid,2)
postmean <- sweep(postmean,2,imod$summary.fixed$mean,"+")
p <- ggplot(reading, aes(cagegrp, piat, group=id)) + geom_line(col=gray(0.95)) + xlab("centered age")
p+geom_abline(data=postmean,intercept=postmean[,1],slope=postmean[,2])

library(gridExtra)
sd.epsilon <- bri.hyper.sd(imod$internal.marginals.hyperpar[[1]],internal=TRUE)
sd.intercept <- bri.hyper.sd(imod$internal.marginals.hyperpar[[2]],internal=TRUE)
sd.slope <- bri.hyper.sd(imod$internal.marginals.hyperpar[[3]],internal=TRUE)
p1 <- ggplot(data.frame(sd.epsilon),aes(x,y))+geom_line()+ggtitle("Epsilon")+xlab("piat")+ylab("density")
p2 <- ggplot(data.frame(sd.intercept),aes(x,y))+geom_line()+ggtitle("Intercept")+xlab("piat")+ylab("density")
p3 <- ggplot(data.frame(sd.slope),aes(x,y))+geom_line()+ggtitle("Slope")+xlab("piat")+ylab("density")
p4 <- ggplot(data.frame(imod$marginals.hyperpar[[4]]),aes(x,y))+geom_line()+ggtitle("Rho")+ylab("density")
grid.arrange(p1,p2,p3,p4,ncol=2)

bri.hyperpar.plot(imod, together=FALSE)

formula <- log(piat) ~ cagegrp + f(numid, model="iid2d", n = 2*nid) + f(slopeid, cagegrp, copy="numid")
imod <- inla(formula, family="gaussian", data=reading)
bri.density.summary(imod$marginals.hyperpar[[4]])
     mean        sd    q0.025      q0.5    q0.975      mode 
 0.059024  0.112174 -0.161646  0.059181  0.276275  0.059951 
data(reading, package="brinla")
reading$id <- as.numeric(reading$id)
newsub <- data.frame(id=90, agegrp = c(6.5,8.5,10.5), piat=c(18, 25, NA))
nreading <- rbind(reading, newsub)
formula <- piat ~ agegrp + f(id, model="iid")
imod <- inla(formula, family="gaussian", data=nreading, control.predictor = list(compute=TRUE))
pm90 <- imod$marginals.fitted.values[[270]]
p1 <- ggplot(data.frame(pm90),aes(x,y))+geom_line()+xlim(c(20,60))
newsub=data.frame(id=90, agegrp = c(6.5,8.5,10.5), piat=c(NA, NA, NA))
nreading <- rbind(reading, newsub)
formula <- piat ~ agegrp + f(id, model="iid")
imodq <- inla(formula, family="gaussian", data=nreading, control.predictor = list(compute=TRUE))
qm90 <- imodq$marginals.fitted.values[[270]]
p1+geom_line(data=data.frame(qm90),aes(x,y),linetype=2)+xlab("PIAT")+ylab("density")

nsamp <- 10000
randprec <- inla.hyperpar.sample(nsamp, imod)[,1]
neweps <- rnorm(nsamp, mean=0, sd=1/sqrt(randprec))
newobs <- inla.rmarginal(nsamp, pm90) + neweps
dens <- density(newobs)
p1 + geom_line(data=data.frame(x=dens$x,y=dens$y),aes(x,y),linetype=2)

Classical Z-matrix Model

library(lme4)
mmod <- lmer(nitrogen ~ 1+(1|site), reeds)
Z <- getME(mmod, "Z")
X <- getME(mmod, "X")
library(INLA)
n <- nrow(reeds)
formula <- y ~ -1 + X +  f(id.z, model="z",  Z=Z)
imodZ <- inla(formula, data = list(y=reeds$nitrogen, id.z = 1:n, X=X))
library(brinla)
bri.hyperpar.summary(imodZ)
                                    mean      sd  q0.025    q0.5  q0.975    mode
SD for the Gaussian observations 0.29079 0.05833 0.19817 0.28283 0.42608 0.26647
SD for id.z                      0.31313 0.15740 0.10893 0.27909 0.71369 0.22175
imodZ$summary.random
$id.z
   ID       mean       sd 0.025quant    0.5quant 0.975quant        mode       kld
1   1 -0.0041086 0.085306  -0.209637  0.00016443   0.156661 -0.00015857 0.0052617
2   2 -0.0041088 0.085306  -0.209637  0.00016426   0.156660 -0.00015874 0.0052617
3   3 -0.0041069 0.085306  -0.209633  0.00016583   0.156664 -0.00015718 0.0052616
4   4 -0.0041101 0.085306  -0.209639  0.00016323   0.156658 -0.00015975 0.0052617
5   5 -0.0041067 0.085306  -0.209633  0.00016596   0.156665 -0.00015705 0.0052616
6   6 -0.0504086 0.147214  -0.530901 -0.00406941   0.025792 -0.00097854 0.0010167
7   7 -0.0504092 0.147213  -0.530901 -0.00407005   0.025791 -0.00098011 0.0010176
8   8 -0.0504093 0.147213  -0.530901 -0.00407012   0.025791 -0.00098028 0.0010177
9   9 -0.0504095 0.147213  -0.530901 -0.00407031   0.025791 -0.00098072 0.0010179
10 10 -0.0504095 0.147213  -0.530901 -0.00407041   0.025791 -0.00098096 0.0010181
11 11  0.0533301 0.154185  -0.027216  0.00339512   0.557732  0.00069156 0.0015510
12 12  0.0533294 0.154185  -0.027216  0.00339500   0.557730  0.00069112 0.0015510
13 13  0.0533319 0.154186  -0.027216  0.00339549   0.557737  0.00069283 0.0015511
14 14  0.0533312 0.154186  -0.027216  0.00339536   0.557735  0.00069237 0.0015511
15 15  0.0533293 0.154184  -0.027216  0.00339496   0.557729  0.00069100 0.0015509
16 16 -0.0041052 0.085314  -0.209677  0.00018201   0.156701 -0.00016665 0.0052413
17 17 -0.0503191 0.147077  -0.530555 -0.00406629   0.025763 -0.00099760 0.0010821
18 18  0.0531454 0.153943  -0.027167  0.00337380   0.557112  0.00069253 0.0016927
data(meatspec, package="brinla")
trainmeat <- meatspec[1:172,]
testmeat <- meatspec[173:215,]
wavelengths <- seq(850, 1050, length=100)
modlm <- lm(fat ~ ., trainmeat)
rmse <- function(x,y) sqrt(mean((x-y)^2))
rmse(predict(modlm,testmeat), testmeat$fat)
[1] 3.814
plot(wavelengths,coef(modlm)[-1], type="l",ylab="LM Coefficients")

n <- nrow(meatspec)
X <- matrix(1,nrow = n, ncol= 1)
Z <- as.matrix(meatspec[,-101])
y <- meatspec$fat
y[173:215] <- NA
scaley <- 100
formula <- y ~ -1 + X +  f(idx.Z, model="z", Z=Z)
zmod <- inla(formula, data = list(y=y/scaley, idx.Z = 1:n, X=X), control.predictor = list(compute=TRUE))
predb <- zmod$summary.fitted.values[173:215,1]*scaley
rmse(predb, testmeat$fat)
[1] 1.9027
rcoef <- zmod$summary.random$idx.Z[216:315,2]
plot(wavelengths, rcoef, type="l", ylab="Ridge Coefficients")

int.fixed <- list(prec = list(initial = log(1.0e-9), fixed=TRUE))
u.prec <- list(prec = list(param = c(1.0e-3, 1.0e-3)))
epsilon.prec <- list(prec = list(param = c(1.0e-3, 1.0e-3)))
idx.X <- c(1, rep(NA,100))
idx.Z <- c(NA, 1:100)
scaley <- 100
formula <- y ~ -1 + f(idx.X,  model="iid", hyper = int.fixed) + f(idx.Z,  model="iid", hyper = u.prec)
amod <- inla(formula, data = list(y=y/scaley, idx.X=idx.X, idx.Z=idx.Z), control.predictor = list(A=cbind(X, Z),compute=TRUE), control.family = list(hyper = epsilon.prec))
predb <- amod$summary.fitted.values[173:215,1]
rmse(predb, testmeat$fat/scaley)*scaley
[1] 1.9013

Generalized Linear Mixed Models

Poisson GLMM

data(nitrofen, package="boot")
head(nitrofen)
  conc brood1 brood2 brood3 total
1    0      3     14     10    27
2    0      5     12     15    32
3    0      6     11     17    34
4    0      6     12     15    33
5    0      6     15     15    36
6    0      5     14     15    34
library(dplyr)
library(tidyr)
lnitrofen <- select(nitrofen, -total) %>% mutate(id=1:nrow(nitrofen)) %>% gather(brood,live,-conc,-id) %>% arrange(id)
lnitrofen$brood <- factor(lnitrofen$brood,labels=1:3)
head(lnitrofen)
  conc id brood live
1    0  1     1    3
2    0  1     2   14
3    0  1     3   10
4    0  2     1    5
5    0  2     2   12
6    0  2     3   15
lnitrofen$jconc <- lnitrofen$conc + rep(c(-10,0,10),50)
ggplot(lnitrofen, aes(x=jconc,y=live, shape=brood)) + geom_point(position = position_jitter(w = 0, h = 0.5)) + xlab("Concentration")

library(lme4)
glmod <- glmer(live ~ I(conc/300)*brood + (1|id), nAGQ=25, family=poisson, data=lnitrofen)
summary(glmod)
Generalized linear mixed model fit by maximum likelihood (Adaptive Gauss-Hermite Quadrature, nAGQ = 25) ['glmerMod']
 Family: poisson  ( log )
Formula: live ~ I(conc/300) * brood + (1 | id)
   Data: lnitrofen

     AIC      BIC   logLik deviance df.resid 
   313.9    335.0   -150.0    299.9      143 

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.208 -0.606 -0.008  0.618  3.565 

Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 0.0911   0.302   
Number of obs: 150, groups:  id, 50

Fixed effects:
                   Estimate Std. Error z value Pr(>|z|)
(Intercept)          1.6386     0.1367   11.99  < 2e-16
I(conc/300)         -0.0437     0.2193   -0.20     0.84
brood2               1.1688     0.1377    8.48  < 2e-16
brood3               1.3512     0.1351   10.00  < 2e-16
I(conc/300):brood2  -1.6730     0.2487   -6.73  1.7e-11
I(conc/300):brood3  -1.8312     0.2451   -7.47  7.9e-14

Correlation of Fixed Effects:
            (Intr) I(c/300) brood2 brood3 I(/300):2
I(conc/300) -0.831                                 
brood2      -0.700  0.576                          
brood3      -0.714  0.587    0.735                 
I(cn/300):2  0.529 -0.609   -0.804 -0.571          
I(cn/300):3  0.538 -0.617   -0.569 -0.804  0.613   
formula <- live ~ I(conc/300)*brood + f(id, model="iid")
imod <- inla(formula, family="poisson", data=lnitrofen)
imod$summary.fixed
                        mean      sd 0.025quant  0.5quant 0.975quant      mode        kld
(Intercept)         1.639578 0.13599    1.36778  1.641091    1.90290  1.644199 6.3039e-12
I(conc/300)        -0.041354 0.21789   -0.47284 -0.040367    0.38417 -0.038459 2.2411e-12
brood2              1.163779 0.13764    0.89730  1.162479    1.43753  1.159878 1.2599e-13
brood3              1.345928 0.13505    1.08493  1.344500    1.61502  1.341638 1.7652e-13
I(conc/300):brood2 -1.663449 0.24847   -2.15557 -1.662005   -1.17976 -1.659093 1.7982e-15
I(conc/300):brood3 -1.820739 0.24499   -2.30621 -1.819229   -1.34406 -1.816188 9.2485e-15
bri.hyperpar.summary(imod)
             mean       sd  q0.025    q0.5  q0.975    mode
SD for id 0.29335 0.057522 0.18809 0.29033 0.41501 0.28449
library(reshape2)
mf <- melt(imod$marginals.fixed)
cf <- spread(mf,Var2,value)
names(cf)[2] <- 'parameter'
ggplot(cf,aes(x=x,y=y))+geom_line()+facet_wrap(~ parameter, scales="free")+geom_vline(xintercept=0)+ylab("density")

bri.fixed.plot(imod)

bri.fixed.plot(imod,together=TRUE)

multeff <- exp(imod$summary.fixed$mean)
names(multeff) <- imod$names.fixed
multeff[-1]
       I(conc/300)             brood2             brood3 I(conc/300):brood2 I(conc/300):brood3 
           0.95949            3.20201            3.84175            0.18948            0.16191 
sden <- data.frame(bri.hyper.sd(imod$marginals.hyperpar[[1]]))
ggplot(sden,aes(x,y)) + geom_line() + ylab("density") + xlab("linear predictor")

bri.hyperpar.plot(imod)

mden <- data.frame(inla.tmarginal(function(x) exp(1/sqrt(x)), imod$marginals.hyperpar[[1]]))
ggplot(mden,aes(x,y)) + geom_line() + ylab("density") + xlab("multiplicative")

formula <- live ~ I(conc/300)*brood + f(id, model="iid")
imod <- inla(formula, family="poisson", data=lnitrofen, control.compute=list(dic=TRUE))
formula <- live ~ log(conc+1)*brood + f(id, model="iid")
imod2 <- inla(formula, family="poisson", data=lnitrofen, control.compute=list(dic=TRUE))
c(imod$dic$dic,  imod2$dic$dic)
[1] 786.11 842.06
mreff <- imod$summary.random$id$mean
qqnorm(mreff)
qqline(mreff)

formula <- live ~ I(conc/300)*brood + f(id, model="iid")
imod <- inla(formula, family="poisson", data=lnitrofen, control.compute=list(cpo=TRUE))
plot(log(imod$cpo$cpo),ylab="log(CPO)")

lnitrofen[which.min(imod$cpo$cpo),]
    conc id brood live jconc
134  310 45     2   10   310
lnitrofen %>% filter(brood == 2, conc==310)
   conc id brood live jconc
1   310 41     2    0   310
2   310 42     2    0   310
3   310 43     2    0   310
4   310 44     2    0   310
5   310 45     2   10   310
6   310 46     2    0   310
7   310 47     2    0   310
8   310 48     2    0   310
9   310 49     2    0   310
10  310 50     2    0   310
pit <- imod$cpo$pit
n <- length(pit)
uniquant <- (1:n)/(n+1)
logit <- function(p) log(p/(1-p))
plot(logit(uniquant), logit(sort(pit)), xlab="uniform quantiles", ylab="Sorted PIT values")
abline(0,1)

which.max(pit)
[1] 134
sdu <- 0.3
pcprior <- list(prec = list(prior="pc.prec", param = c(3*sdu,0.01)))
formula = live ~ I(conc/300)*brood + f(id, model="iid", hyper = pcprior)
imod2 = inla(formula, family="poisson", data=lnitrofen)
bri.hyperpar.summary(imod2)
             mean       sd  q0.025    q0.5  q0.975    mode
SD for id 0.30966 0.056721 0.20718 0.30634 0.43046 0.30016
lnitrofen$obsid = 1:nrow(nitrofen)
formula <- live ~ I(conc/300)*brood + f(id, model="iid") + f(obsid, model="iid")
imodo <- inla(formula, family="poisson", data=lnitrofen)
bri.hyperpar.summary(imodo)
                 mean        sd    q0.025      q0.5   q0.975     mode
SD for id    0.293850 0.0564034 0.1951276 0.2899458 0.415677 0.283363
SD for obsid 0.010303 0.0062123 0.0037649 0.0085252 0.027042 0.006218

Binary GLMM

data(ohio, package="brinla")
table(ohio$smoke)/4

  0   1 
350 187 
xtabs(resp ~ smoke + age, ohio)/c(350,187)
     age
smoke      -2      -1       0       1
    0 0.16000 0.14857 0.14286 0.10571
    1 0.16578 0.20856 0.18717 0.13904
library(lme4)
modagh <- glmer(resp ~ age + smoke + (1|id), nAGQ=25, family=binomial, data=ohio)
summary(modagh)
Generalized linear mixed model fit by maximum likelihood (Adaptive Gauss-Hermite Quadrature, nAGQ = 25) ['glmerMod']
 Family: binomial  ( logit )
Formula: resp ~ age + smoke + (1 | id)
   Data: ohio

     AIC      BIC   logLik deviance df.resid 
  1603.3   1626.0   -797.6   1595.3     2144 

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.373 -0.201 -0.177 -0.149  2.508 

Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 4.69     2.16    
Number of obs: 2148, groups:  id, 537

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -3.1015     0.2191  -14.16   <2e-16
age          -0.1756     0.0677   -2.60   0.0095
smoke         0.3986     0.2731    1.46   0.1444

Correlation of Fixed Effects:
      (Intr) age   
age    0.244       
smoke -0.493 -0.008
formula <- resp ~ age + smoke + f(id, model="iid")
imod <- inla(formula, family="binomial", data=ohio)
imod$summary.fixed
                mean      sd 0.025quant 0.5quant 0.975quant     mode        kld
(Intercept) -2.99276 0.20434  -3.416832 -2.98475  -2.614898 -2.96902 2.1272e-14
age         -0.16665 0.06283  -0.290606 -0.16646  -0.043901 -0.16606 1.8456e-14
smoke        0.39136 0.23954  -0.077759  0.39068   0.863758  0.38935 1.3287e-12
ilogit <- function(x) exp(x)/(1 + exp(x))
ilogit(imod$summary.fixed[1,c(3,4,5)])
            0.025quant 0.5quant 0.975quant
(Intercept)   0.031774  0.04812   0.068186
exp(imod$summary.fixed[-1, c(3,4,5)])
      0.025quant 0.5quant 0.975quant
age      0.74781  0.84666    0.95705
smoke    0.92519  1.47798    2.37206
bri.hyperpar.summary(imod)
            mean      sd q0.025   q0.5 q0.975   mode
SD for id 1.9258 0.16156 1.6253 1.9196 2.2603 1.9097
exp(bri.hyperpar.summary(imod)[3:5])
[1] 5.0800 6.8180 9.5858
table(xtabs(resp ~ id, ohio))

  0   1   2   3   4 
355  97  44  23  18 
library(gridExtra)
p1 <- ggplot(data.frame(imod$marginals.fixed[[1]]),aes(x,y))+geom_line()+xlab("logit")+ylab("density")+ggtitle("Intercept")
p2 <- ggplot(data.frame(imod$marginals.fixed[[2]]),aes(x,y))+geom_line()+xlab("logit")+ylab("density")+ggtitle("age")
p3 <- ggplot(data.frame(imod$marginals.fixed[[3]]),aes(x,y))+geom_line()+xlab("logit")+ylab("density")+ggtitle("smoke")
sden <- data.frame(bri.hyper.sd(imod$marginals.hyperpar[[1]]))
p4 <- ggplot(sden,aes(x,y)) + geom_line() + xlab("logit") + ylab("density")+ggtitle("SD(u)")
grid.arrange(p1,p2,p3,p4,ncol=2)

inla.pmarginal(0,imod$marginals.fixed$smoke)
[1] 0.050916
formula <- resp ~ age + smoke + f(id, model="iid")
imod <- inla(formula, family="binomial", data=ohio, control.compute=list(dic=TRUE,waic=TRUE))
formula <- resp ~ factor(age) + smoke + f(id, model="iid")
imod1 <- inla(formula, family="binomial", data=ohio, control.compute=list(dic=TRUE,waic=TRUE))
formula <- resp ~ factor(age)*smoke + f(id, model="iid")
imod2 <- inla(formula, family="binomial", data=ohio, control.compute=list(dic=TRUE,waic=TRUE))
c(imod$dic$dic, imod1$dic$dic, imod2$dic$dic)
[1] 1466.3 1463.2 1463.5
c(imod$waic$waic, imod1$waic$waic, imod2$waic$waic)
[1] 1418.9 1416.3 1417.7
exp(imod1$summary.fixed[-1, c(3,4,5)])
              0.025quant 0.5quant 0.975quant
factor(age)-1    0.74517  1.08375    1.57718
factor(age)0     0.65627  0.95965    1.40216
factor(age)1     0.38379  0.57802    0.86375
smoke            0.92488  1.48190    2.38575
ohio$obst = rep(1:4,537)
ohio$repl = ohio$id + 1
formula <- resp ~ factor(age) + smoke + f(obst, model="ar1", replicate = repl)
imod <- inla(formula, family="binomial", data=ohio)
exp(imod$summary.fixed[-1, c(3,4,5)])
              0.025quant 0.5quant 0.975quant
factor(age)-1    0.74421  1.08461    1.58185
factor(age)0     0.65311  0.95868    1.40595
factor(age)1     0.37896  0.57416    0.86227
smoke            0.92547  1.48607    2.39799
bri.hyperpar.summary(imod)
                mean       sd  q0.025    q0.5  q0.975   mode
SD for obst  1.97258 0.170482 1.66238 1.96343 2.33132 1.9443
Rho for obst 0.98471 0.014733 0.94575 0.98887 0.99856 0.9966
cmod <- inla(formula, family="binomial", data=ohio, control.inla = list(correct = TRUE))
cmod$summary.fixed
                   mean      sd 0.025quant  0.5quant 0.975quant      mode        kld
(Intercept)   -3.114881 0.28248  -3.708581 -3.101354   -2.59882 -3.075079 2.7827e-14
factor(age)-1  0.088258 0.19950  -0.302904  0.088034    0.48030  0.087609 1.5160e-13
factor(age)0  -0.045929 0.20463  -0.448273 -0.045802    0.35532 -0.045534 5.5515e-13
factor(age)1  -0.600099 0.22165  -1.041175 -0.598184   -0.17010 -0.594445 5.7201e-13
smoke          0.430558 0.26989  -0.096802  0.429313    0.96421  0.426899 1.7247e-12
bri.hyperpar.summary(cmod)
                mean       sd  q0.025    q0.5  q0.975    mode
SD for obst  2.33088 0.218061 1.90502 2.33241 2.75712 2.34995
Rho for obst 0.96272 0.024311 0.90517 0.96735 0.99461 0.98362
sessionInfo()
R version 3.4.0 (2017-04-21)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.5

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] reshape2_1.4.2  tidyr_0.6.1     dplyr_0.5.0     gridExtra_2.2.1 ggplot2_2.2.1   brinla_0.1.0    INLA_17.06.30  
 [8] sp_1.2-4        lme4_1.1-13     Matrix_1.2-9    rmarkdown_1.4  

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.10     compiler_3.4.0   nloptr_1.0.4     plyr_1.8.4       tools_3.4.0      digest_0.6.12   
 [7] evaluate_0.10    tibble_1.3.0     nlme_3.1-131     gtable_0.2.0     lattice_0.20-35  DBI_0.6-1       
[13] yaml_2.1.14      stringr_1.2.0    knitr_1.15.1     rprojroot_1.2    grid_3.4.0       R6_2.2.0        
[19] minqa_1.2.4      magrittr_1.5     backports_1.0.5  scales_0.4.1     htmltools_0.3.5  MASS_7.3-47     
[25] splines_3.4.0    assertthat_0.2.0 colorspace_1.3-2 labeling_0.3     stringi_1.1.5    lazyeval_0.2.0  
[31] munsell_0.4.3